Equation of states for classical Coulomb systems. 
The use of Hubbard-Schofield approach. 

J. Ortner 

Institut fur Physik, Humboldt Universitat Berlin 

Invalidenstr. 110, 10115 Berlin, Germany 

(to be published in Phys. Rev. E) 

V»L', An effective method based on Hubbard-Schofield approach [Phys. Lett. A 40, 245 (1972)] is 

ji. . developed to calculate the free energy of classical Coulomb systems. This method significantly 

fp. ' simplifies the derivation of the cluster expansion. A diagrammatic representation of the cluster 

integrals is proposed. Simple rules providing the leading order in density n of each diagrammatic 
contribution are found. We calculate the n 3 contribution and recover the results at the order n 5 ' 2 
obtained by the traditional method of resummation of diverging Mayer bonds. 
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H, 1 I. INTRODUCTION 

This paper is adressed to the study of the virial expansion of the Helmholtz free energy (hereafter the free energy) 
r ■ for a classical Coulomb system. We consider a multicomponent system of point-like ions embedded in a neutraliz- 
ing background. Because of the long-range Coulombic nature of the interaction potential between two charges the 
corresponding virial expansion Mayer graphs pj diverge. In the traditional method these long-range divergencies are 
removed via the chain resummations first introduced by Mayer J2J and Salpeter [|3| and further developed by the works 
of Meeron B and Abe pj. To avoid complicated calculations related to the chain resummation we present in this 
^"v paper an alternative method of calculation of classical Coulomb systems thermodynamic functions. 

We show in this paper that the earlier results can be obtained much easier using the method of collective variables 
and integral transformations. Moreover, this method makes it possible to obtain the free energy of a classical plasma 
system in a systematic manner up to an arbitrary order (at least in principal). In addition, this method is capable 
to describe not only the low density limiting thermodynamic behavior of a Coulomb system but also the region of 
critical density of Coulomb fluids and that of a strongly coupled plasma. However, in this paper we will restrict our 
considerations to the low density region of a classical ion mixture. 

The method of collective variables is a powerfool tool of investigation of both classical and quantum Coulomb 

systems. There were two basic lines in this method mainly developed already in the fifties. The first line starts with 

f~^ ■ the initial plasma Hamiltonian being the sum of charges kinetic energy and the interaction potential between them 

0^ ' and converts it into the collective variables Hamiltonian using the canonical transformations H . 

0^ . The second line j7||§ starts from the configurational integral of individual particles. Certain transformations lead 

r/3 ' then to the configurational integral expressed through collective variables. 

Using the first approach the physical processes can be interpreted in terms of the collective variables. Thus the 

plasma oscillations are represented as oscillations of the density fluctuations Fourier components. Following the second 

^~>j general line it is possible to obtain an expansion of the cluster integral in a systematic manner. 

i-^ Using the formal representation of the configurational integral in terms of the collective variables Zubarev J7J and 

J~7^ Juchnovskij gave a cluster expansion of the free energy. The transformation from the individual (i.e., the position of 

£> . the particles) to the collective coordinates (i.e., the Fourier transforms of the density fluctuations) is performed via the 

corresponding Jacobian. Within the Random Phase Approximation for the Jacobian one arrives at the Debye-Hiickel 

approximation for the free energy J7J , and a systematic expansion of the Jacobian leads to the cluster expansion of 

C^ ■ the free energy f|fl. However, this method is still cumbersome. 

By use of integral transformations the method of collective variables have become technically more feasible. Here, 
instead of a Jacobian transformation an identity is used which expresses the Coulomb interaction in terms of external 
interactions pQ-Q. As for Coulomb systems this exact version of the mean field idea of Debye-Hiickel leads to 
the sine-Gordon (SG) representation of the configurational integral p^WLl. By virtue of a rigorous mathematical 
pro of it was shown that the SG transformation produces a cluster expansion of a system with long-range interactions 
ft5| . Further the SG theory of Coulomb gas was used to analyze the metal-insulator transition |U|[L7) and Coulomb 
criticality Jl8f . Such a way of analysis of the Coulomb criticality was claimed to be the most promissing [pi . However, 
the SG representation in the pure Coulomb version is applicable only for lattice models (with possibly vanishing lattice 
constant) and for pointlike charges. 



A hybrid method combining the advantages of the Zubarev-Juchnovskij approach with that of Stratonovich, Hub- 
bard and Kac was developed in Ref. P0(| . Here the particle interactions are divided into long- and short-range 
interactions. Thus short-range repulsive interactions as, e.g., hard core repulsions are introduced in a natural way. 
Using the Stratonovich- Hubbard-Kac transformation the original system with the both types of interactions is mapped 
onto a reference system with short range interactions only. In contrast to the SG theory the Hubbard-Schoficld (HS) 
approach is capable to perform off-lattice calculations. Further in the HS approach the cumulant expansion is used 
to map the Hamiltonian of the original system (a nonionic or ionic fluid) onto an effective magnetic-like Hamiltonian. 
Since the magnetic system Hamiltonian can be recasted into the Landau-Ginzburg- Wilson form |23[| such a represen- 
tation is very convenient for the analysis of both ordinary fluids |E1[ and charged hard sphere systems criticality |E2) . 
To our knowledge in Ref. |g4| the HS approach was first applied to a one-component plasma. 

Despite the significant success of the method of collective coordinates in describing Coulomb systems the low 
density limit of Coulomb systems was mainly studied by traditional methods of statistical physics p5[ . Starting by 
the charging formula for the free energy an expansion in terms of ordinary Mayer functions can be found. However, in 
the case of Coulomb systems the Mayer functions diverge. Expanding the Mayer function in powers of the Coulomb 
potential one can collect the Mayer series into special subseries and perform a partial summation of Mayer bonds. 
The resummed Mayer bond is representable through the screened Coulomb potential and is integrable. Thus the sum 
of all ring diagrams gives the Debye-Hiickel approximation for the free energy. Picking up in all diagrams the chain 
of Mayer bonds leading to the screened potential and ordering all diagrams by the number of vortices the cluster 
expansion for the free energy is obtained. The virial expansion up to the order n 5 ' 2 (n being the number density 
of ions) was found within this method k\. 20j27],E5| . However, to perform this procedure a number of very refined 
diagrammatic transformations are required ||25|. Various first-principles formalisms based on the chain resummation 
of Mayer bonds have been used to generalize the virial expansion of classical Coulomb systems to the case of quantum 
plasmas. First, the method of effective potentials introduced by Morita |28| for quantum systems, has been applied by 
Ebeling to the Coulomb case |29| . A renewed interest in the exact calculation of thermodynamic functions for weakly 
coupled quantum plasmas beyond the Debye-Hiickel limiting law has been observed recently. The virial expansion 
for the free energy up to the order n 5 ' 2 in the density was derived using the Feynman-Kac formalism ||30| and the 
method of Green's functions |3l{ ], respectively. The latter method involves an additional expansion with respect to 
the square of charge e 2 . Confirming and completing the earlier results by Haga Pq| , Friedman Pq] , DcWitt |32|] and 
Ebeling E9| the both approaches arrived at the same result. Therefore the virial equation of state for both classical 
and quantum Coulomb systems is now generally accepted as known up to order n 5 ' 2 . 

In this paper we show that the density expansion of a classical Coulomb system can be obtained much easier 
using the method of collective coordinates. As Zubarev and Juchnovskij we rewrite the partition function in terms 
of the collective variables. In contrast to them we will not use the Jacobian transforming the space coordinates of 
individual particles into the collective variables. Instead, we represent the configurational integral of the Coulomb 
system through the configurational integral of a noninteracting system in external fields. Further following the line of 
Hubbard and Schoficld we employ the cumulant expansion. This enables us to rewrite the configurational integral in 
a similar way as the partition function for the magnetic system having an Ising-like Hamiltonian. The coefficients of 
the effective magnetic-like Hamiltonian are expressable in terms of the ideal gas correlation functions and are perfectly 
known. Performing a perturbation expansion in terms of the anharmonic contributions to the effective Hamiltonian 
one directly obtains the cluster expansion of the configurational integral and of the free energy of a classical plasma 
without bonds resummation (Section |l|). Instead the bare Coulomb interaction is automatically screened in the HS 
approach. Ultimately, the integrals that determine the virial coefficients in the HS method are the same which appear 
in the more traditional Abe-Meeron-Friedmann approach. In this sense the HS method may be considered as another 
evidence of the fact that all virial coefficient are expressable in terms of Mayer-like graphs built with screened bonds. 
Arriving at the virial equation of state it will be shown that the method of collective variables in the Hubbard-Schofield 
representation is capable to describe not only the critical behavior but also the low-density limit of Coulomb systems. 

Additionally, there are some basic ideas in the literature on how to generalize the method of collective variables to 
the case of a strongly coupled plasma. We mention here the calculation of "first-principles" expressions for the free 
energy over the entire density region of classical plasmas [Bq,M ■ In this approach the static structure factor in the 
Debye-Hiickel form is employed, strong coupling effects are involved by introducing an upper bound for the collective 
variables wavevector, as it was done by Debye in the theory of specific heats of solids. 

Further approaches are devoted to the study of dynamic properties of coupled Coulomb systems using the repre- 
sentation of the plasma Hamiltonian into collective variables. In a rather incomplete list we mention the approaches 
in Refs. |34||35| based on the Mori-Zwanzig theory |37],|3(| (or the memory function formalism), the approaches based 
on the theory of moments | B8J , and the approach based on the quasilocalized charge approximation 39 1 . 

The present paper is organized as follows. In Sec. |J, we a PPly t ne approach of Hubbard and Schoficld |2(]] to 
map the original Hamiltonian of the classical Coulomb system onto a magnetic-like Hamiltonian. The coefficients 
of the obtained magnetic Hamiltonian are expressed via the ideal gas structure factors, and are calculated exactly. 



Expanding the anharmonic contributions of the magnetic Hamiltonian we obtain the cluster expansion of a classical 
plasma. 

The cluster expansion is the starting point for the density expansion of the classical Coulomb system. Explicit 



calculations are carried out up to the order n in Sec. [II. Comparison to previous results is carried out briefly 



II. CONFIGURATIONAL INTEGRAL 

Consider a classical plasma consisting of M sorts of ions with 7V ions of a given species a with masses m a and 
charges e a . In what follows we will omit the summation bounds in the summation over the particle types a if the 
summation is carried out from 1 to M . The total number of particles is N — ^2 a N a . The system of ions is immersed 
into a neutralizing background. The plasma system with total volume V and temperature in energy units T — 1//3 is 
described by the interaction potential 



u{fi,...,fN) = ^EEE^ 6 ^' 7 *^ ■ w 

a.b i—lj — 1 

As far as we consider point-like ions the interaction potential between two particles is given by the Coulomb potential 

v ab (n,fj) = p^ . (2) 

Introduce now the collective variables of the plasma system, the Fourier transforms of the charge density, 

en = (-r 1/2 E e - E ^ ■ ( 3 ) 

where v = Y, a e 2 a N a , and N a = Ya,=i N b- 

The interaction potential can be expressed in terms of the collective variables, 

PU=lj2a(k)[g j :g_ j : - l] , (4) 

fc#0 

with the notation a(k) — n 2 /k 2 and k 2 = (4w)/(VT) J^ a e a-^a being the square of the inverse Debye radius. The 
contribution with k — cancels due to the presence of the background. 
Our aim now is to calculate the free energy of the plasma system, 

F = - T In Z ld Z c . (5) 

Here Z^ is the ideal part of the partition integral, whereas the configurational integral is given by 

Z c = ^ Jdfi...drZe-l 3u <*-> r V . (6) 

The configurational integral is expressed now through the collective variables (?]|| 



ex P I 2 E a{ ^> ( exp I ~ 2 51 a ( fc ) g k g -k 

\ fe#0 / \ \ fc#0 / 



(7) 



where the angular brackets denote averaging with respect to an ideal gas system, i.e., 



(...) = ^L /^...,/m (...) . (SI 



Using the identity 



i„v) . (wr* r^(-\i 



exp(--a 2 a; 2 ) = (27ra 2 ) '" \ cxp ( -- y — ) exp ( - ixy) dy (9) 



one expresses the configurational integral of the interacting system through the configurational integral of a system 
of noninteracting particles moving in an external field (p^ |2C|] 



z c = CX pQ]>>(fc))n{[ 2 ™(fc)r 1/2 /<%} 



exp 



-(Va(fc) 1 <p#p_ i 



l J2k e sv>-k 



(10) 



Note that the field variable ipj: (as well as g^) is a complex number, tpj: — (^g c + i if^ s with real <pg c and tpj: s - Therefore 
in Eq.(|lO|) the product of integrals has to be understood in the following manner: 



dtpzdcpz 



OO /"DC 



d( Pkc d( Pks and <P#P-i = Wkc 



\Vk s 



— oo J — OO 



In what follows we will use the shorter notations of Eq. ( [i(j|) 
Regard the functional 



*M = (cxp< 1^2,6% 



<P-% 



The moments of g^ (i.e., the structure factors of an ideal gas) can be obtained by functional differentiation 

d s t>[tp] 



Consider further the functional 



H) 6 



d<p_ %i d<p_ %2 ...dip_ le 



tp=0 



Q[y] = ln*[ 
The Taylor expansion of defines the cumulants of g [[fOJ 

e 



9 M = 5Z "7 5Z u n(kl, fo, ...,k n 



<P-% 1 <P-% 2 ---'P-% n 



k\ ...k n 



i.e., the cumulants of order s can be obtained by the functional derivative 

<9 S 6M 



u s (ki, k 2 , .. ., k 3 ) = (-i) 6 



&P-lf<P-%>---&P-k. 



<p=0 



The first four cumulants read 

ui (fc) = (q s ) =0 , k / , 



«2 



(11) 



(12) 



(13) 



(14) 



(15) 



(16) 



(fci,fc 2 J = (^0*5) > fci ^ o,fc 2 ^ o 

U3(fci,fc2,fc 3 J = i^g^g^g^J , fci 7^ 0, fc 2 7^ 0, fc 3 ^ , 

- \Q^Q^) (Ok^k,) ~ (e&8&) (^^3)' fci ^ 0, ... , fcl 7^ . (17) 

In general one expresses the cumulants of order s u s I fei, Ak, ■ • • , k n ) through the ideal gas structure factor of order 
s ( g fe - g fe - ■■• (?£* ) minus all combinations of structure factors of order less than s. 
Calculate now the ideal gas structure factor of order n {k\ 7^ , . . . , k n 7^ 0) pl|. 



e*i«jfi •••«*; 



Na ± +N al M an +N a 



If * \ * " " 

v ~ n,2 yN I df>1 ■■■drN Yl e «l e «2 ■ ■ ■ e a„ 51 ■ ■ ■ X 



1=J V ai +l i„=A/" a „+l 



exp 



i i (fcir^ 



rC^i / 7 " 



0} 



In the thermodynamic limit (i.e., neglecting terms of the order 1/iV) we obtain for the cumulants 



(ki, k 2 , ..., k„) 



u n (fei, fc 2 , . . . , k n ) = N 2 c n 5,- 



ki+...+k„ ,0 



X)** ^ 0» 0C{l,2,...,n} 



ieQ 



with (5 a (, being the Kronccker delta, and 



{N/v 



>-2)/2 , -1 



E e "^ 



(18) 



(19) 



(20) 



The conditions in Eq.(U9) reflect the facts that the sum of all n wavevectors being the argument of 

u n ( k\, fc 2 , . . . , k n \ equals zero, whereas all the sums of wavevectors containing less than n wavevectors k t have 

not to be equal to zero. 

Inserting Eq. (fi9|) into Eq. (p"5|) , and using the definition of the functionals 8 (Eq. dl3) and V& (Eq. (JT3) , we obtain 



/ 00 .^ 

^-^ \ I it"' / (n — 2) 

^ n euip ~n = cx p X, v\ zZ N 2 c « *£+...+£ , v-up-% ■ ■ ■ ?_*, 



(21) 



k\ ...k n 



where the prime at the summation sign indicates the truncation, i.e., the wavevectors ki are choosen in such a way 
that all the sums of less than n wavevectors have not to be equal to zero. 

Substituting Eq.rt2^) into Eq.([l0|) we obtain the Hubbard-Schofield representation of the configurational integral 
for the classical ion mixture 1 20 



Z c = cxp(-Y,a(k))l[\l2na(k)} 
with the effective Hamiltonian 



-1/2 



dtpA e 



■H[<P] 



(22) 



and 



HVp\ = E H w n (ki, k 2 , ..., k n ) ?-%!?-%,■•>?-%„ > 



n=l C, 



k± , k n 



Wi = 



W 2 = 2 S k 1+ k 2 ,0 



«C 3 r 

w 3 — r o,r 



h 



3! N 



I u k 1 +k 2 + k 3 ,0 ' 



(23) 



% Cn 



"V, 



6c 



n'.N 2 



We rewrite the configurational integral in the form 

' exp (a(fc)/2) 



z c = U 



s/2ira{k) 



dipt exp 



--(a \k) + \)<p#>_ % 






(24) 



(25) 



where the 7i„ are defined by, 



TLr, 



i\N~ 



E ' ** 



+fe 2 +...+fe„ ,o 



V-ifc, <P-% 



k 2 



•V-L. 



(26) 



ki_...k n 



with c ra from Eq.(|20|). The Ti n contain the anharmonic contributions of the auxiliary field variables p^ into the 
configurational integral. Neglecting the anharmonic contributions one arrives at the Debye-Huckel approximation for 
the free energy first derived within the framework of collective variables by Zubarev (7) (see the next section) . Notice 
that the presence of the Debye potential in the covariance of the Gaussian measure over the </>g's is the reason for the 
automatic screening within the present approach. 

In what follows in order to obtain the cluster expansion within the HS approach we have to perform manipulations 
analogous to those in the Zubarev- Juchnovskij approach. Expanding the exponential exp {— X)^°=3 H-n} one represents 
the configurational integral as 



jj «p(a(*)/2) J 


-~(a \k) + 1)<P0>_ & 


■ {1 + C[ + C 2 + C 3 + C' A + . . .} , 


abbreviations 


C 1 = H 3 + Ha + 7i 5 + . . . , 


c 'i = \ {[^3 n 3 + Hi Hi + ...] + 2 [n 3 Hi + n 3 n 5 + ... + n± 


(-3 = ^{[^3^3% + H4K4W4 + ■••] + 3 [H 


]TL 3 TLa + TL 3 TL 3 H5 + 



are introduced. Eq. 



TiiTL^TL^ + . . .] +6 [H 3 Tii TL^ + Ti 3 Ti^ 7i 6 . . .]} 
is an infinite sum of products of integrals of the type 



dipt exp 



-\{a \k) + l)<p#p_ n 



[<Pi\ 



f-i 



(27) 



]} 



(28) 



(29) 



The result of the integration over p-r differs from zero only if n = m in Eq.(p9|). Having this in mind and taking into 

account that in TL n the sum over the fe's is truncated, one finds that all linear powers of TL n (i.e., C[) and also a lot of 
products of Ti n hi the higher order C' n (n > 1) vanish under the integration over ip?. In general we have to choose 
products of 



w. 



ll I "'11 j • ■ ■ 1 ""lm J • • • W ng I K s \ , k srls J 



V-fc^-E. 



-V>-%, 



<p-%., <P-z. 



■<p-%. 



(30) 



(see Eqs.( [23| , f24| )) such that one can pair all p-^ with other <p^ whereby i ^ I. 

Consider as an example the pairing of field variables in the lowest nonvanishing product (1/2)7^3 7^3 contributing 
to C' 2 - In this case the integrand is 



n 

k 

2 3! 



exp (a(fc)/2) 

^2?ra(fc) 



dip% exp 



-(a \k) + l)<p £ ip_ U 



n 



jc 3 n 1 1 -c 3 n i J2 (n^n, 

ki,k 2 ,k 3 

cxp(a(fc)/2) 1 (-1) 3 



¥>_£, V-to ¥>-2, 



^1+^2+^3,0 



yMW 



1 



2-3! 



N~'& 



fel , &2 ;ks 



"kl_+%2+k*3 ,0 



-Hfci) + 1] [a- 1 (k 2 ) + 1] [a-Hto) + 1] 



(31) 



It was taken into account here that one has 3! possibilities to pair the field variables p^ in the product 7i 3 Ti 3 . Notice 
also that in the r.h.s. of Eq. (J3l|) the truncation in the sum over k\ , k 2 , k 3 is omitted. Performing in Eq.d27n the 
integration over the field variables ip* one arrives at the following formula for the configurational integral 



*-n{ 



exp (a(fe)/2) 

V a ( k ) + 1 



{1 + Ci + C 2 + C 3 + C 4 



•} 



(32) 



where the C n (without prime) denote the sum of all products of n functions Hi (i = 3 , 4 , . . .) after integration over 
the field variables. Having in mind the above rules for the pairing of field variables under the integration sign one 
finds that all the contributions into C\ vanish, the second cluster C2 is composed of the terms (H n H n ), n > 3 only, 
whereas in the third cluster C3 only the combinations (H n H m Hi) (3 < n < m < I) with n + m > I and n + ra + I 
being an even number remain. We obtain the following structure for the first three clusters 

Ci = , 

C2 = ^{{(H 3 H 3 ) + (H 4 H 4 ) + ...]} , 

C3 = 7^ {[{H4H4H4) + (HqHqHq) + . . .] +3 [(H3H3H4) + [HaHaHq) + 

+ {HiH 5 H 5 ) + ...] +6 [(HsHiHs) + (H 3 H 5 H 6 ) . . .]} . (33) 

Here (. . .) means that the integration over the field variables ip was performed. 

Following Juchnowskij [91 we introduce a diagrammatic representation of the expressions appearing in the clusters 
C n - The diagrams consist of lines and points. The connection of lines and points is called a vortic. An unpaired H n 
contribution occuring in the clusters C' n can be symbolized by a vortic with n open lines 

1 

2 



11 

whereas a paired H n H n contributing to the clusters C n can be represented by the following diagram, 

1 




n 

Comparing the diagrammatic representation with the analytic expression from Eq.(pl|) one can establish the fol- 
lowing rules for the diagrams: 

1 . Outgoing lines are associated with a wavevector ki , incoming with — k, t . 

2. A vortic with outgoing lines 

n 

represents the factor 

ATI 

1 ^n c 

0,7 



N 2 



3. n lines connecting two points replace the notations 

1 



n\ [a- l {k{] + 1] [a- l {ki) + 1] ... [^(M + 1] ' 

4. One has to count the number of topological identical diagrams which differ only by the numeration of the points. 
(In what follows this number will be shown explicitely.) 

5. Summation over the wavevectors has to be performed. 

With the thus defined diagrams one easily expresses the clusters C n . From Eq.(|33]) and calculations analogous to 
that in Eq.(|31[) one finds that the cluster C n is given by the sum of all possible diagrams satisfying the following 
conditions: 

1. The diagrams consist of n points and closed lines. 

2. The minimal number of outgoing and incoming lines in each point of a corresponding diagram equals three. 

3. Selftruncation is forbidden, i.e., closed loops are absent, each line connects two different points of a diagram. In 
addition a connected diagram or a connected subdiagram of a disconnected diagram should have the following 
property. By cutting an arbitrary number of connections of a certain (but arbitrary) point with other points 



but keeping at least one connection one remains with a connected diagram (or subdiagram) still. (This rule 
represents the prime in Eq.([26|).) For example due to this rule the following diagrams should be excluded: 




Since our final aim is the calculation of the free energy we rewrite Eq.(32) in a more convenient form, 

S 3 + Si + ...} , (34) 



Zc = n |exp( # V|| exp{52 



y/aJk) 



1 



The cluster S n can be expressed by the same diagrams as the cluster C n with one additional condition: 

4. From all the diagrams contributing to S n we have to choose only irreducible diagrams, i.e., those which do not 
represent a combination of lower order diagrams. (The latter arise in Eq.(p2) by expanding the exponent in 
Eq.(p4[)). For example the diagram 



contributes to d but does not contribute to 1S4. 
The first cluster integrals S n will be considered in the next section. 



III. VIRIAL EXPANSION 

Our aim now is to obtain the virial expansion for the free energy of the ion mixture. The basic expression that is 
used as a starting-point is the cluster expansion for the interaction part F c of the Hclmholtz free energy per volume 
V in units of T = 1/(3, 



where 



Tdh 



PF C _ 



vt*« 



V 



±J2 W fc ) -- ln ^) + 1]} 



K 

12tT 



(35) 



(36) 



is the leading Debye-Huckel contribution to the free energy M. 

In the following part of this section we will consider the first cluster integrals and extract the leading terms of the 
density expansion within each cluster integral. 



A. Second cluster integral 

The second cluster integral is given by the sums of diagrams 

1 , 
S 2 

or in analytic representation 

S 2 



1 ^ (-l)"c% f dfa dk 2 ... dk n 

9. 2.^ n\N n -' 2 



2 ^ n 

n— 3 



(2tt) 3n 





h 


!+k 2 +... 


ffcn 


.0 




' k 2 
^ + 1 

K 2 




' k 2 

% + 1 

K 2 






^ + 1 



(37) 



(38) 



Substituting the c„ values from Eq. (EG) , and using the convolution theorem and the expression for the inverse Fourier 
transformation of 



d k ^irl a 



b ikr 



I \ '•ab - h 

g ab (r) = — e 
r 



(2nfk 2 + k 
we find the familiar expression for the second cluster integral 



Lb = e a e b P 



(39) 



* 2 = £^E^/<^. ■"«<■" >r 



= E 



2V 

a,b n— 3 

N a N b 



a . b 



2V 



dr 



®ab{r) ~ ^9 2 ab (r) 



where the generalized Mayer function 



$„ 



e~ 9ah{r) - 1 + g ab (r) , 



(40) 



(41) 



is introduced. 

Thus performing the inverse Fourier transform one can establish the following rules of evaluation of a cluster 
expansion graph in the coordinate representation: 

1. Each point is associated with a particle of species a, and provides an integration over the coordinate space - 
nc a J dr a , the total particle density n — N/V and the particle concentration of species a, c a = N a /N, are 
introduced here. 

2. I solid lines connecting the points a and b give a factor ((— 1) //!) [gab(r a b)) 1 where r ab = \r a — r b \. 

3. Introducing a graphic representation of the generalized Mayer function, 



CD 



(42) 



one associates a bond 

exp[- g a b(r a b)} - 1 + g a b(r a b)- 



connecting the points a and b with the Mayer function <frab(r a b) 



4. Finally, summation over the particle types has to be performed. 
Given this rules one represents the second cluster integral in the following manner 



(43) 



In order to get the density expansion of the second cluster integral we expand the corresponding expression in 
Eq.([40|) in powers of k. Up to the order n 3 we get, 



S 2= 2nJ2^ h'rr 2 



a , b 



v 



i( l f ] (r ; " 



1) + Mr) 



U l -f) (, •• ■ I +4m-) +/„,,/.•/„,,(» 



(44) 



' ( '"'' Y ( - 5 - - 1 + 5«r - ^-nr) + Ul ab nf f ab (r) - l l ab n\ ( f ab (r) + \( 1 -^ " 



5! \ r 
where the truncated bond 



6 V t 



r , ■ I a-b \ -. , Lab *■ I ^ab 

/a6 (r) = exp(--j -l + _--(_ 



(45) 



has been introduced. 

Similar calculations of the second cluster integral at the order n 3 for a quantum plasma have been performed by 
Kahlbaum |42) and Ortner Q. 

A straightforward calculation of the integrals in Eq.(B3) leads to the following result, 



S 2 = vL^CaCbll 
^ a,b 

+ n 2 K 2 y^c a c b l 5 a 



a , b 



Innlab + a-2 



57r n , 

— Intvlab + C 2 



Tl K y ^ CaCbtab 

a . b 



111 Klab 



(46) 



with the numerical constants 



"2 



C-2 



= 3 


2C B 


6 


= 0.439519 


7T 

3 


2C £ - 


6 


= - 0.30641 


57T 

24 


2C B - 


203 , . 

hln5 

60 


= - 0.4054 



(47) 



Ce being the Euler-Mascheroni constant, Ce = 0.57721566 . . .. 

Thus we recover the result of Abe p] (i.e., the density expansion of the free energy up to the order n 2 ) and have 
written down all terms of the order n and n 3 In n stemming from the second cluster integral. 



B. Third cluster integral 



Following the general rules of construction of the cluster integrals of arbitrary order one expresses the third cluster 
integral through the following sum of diagrams M 



*-5i ' 







(48) 



With the graphic representation of the generalized Mayer function we represent the third cluster integral in the 
compact manner fBPqM4], 



& 




+ 




(49) 



Using the diagrammatic rules one easily writes down the explicit expression for the third cluster integral 
,3 



Sz = —^^c a CbC c ldr a dr b dr c - 3g afc (r Q6 )$ QC (r QC )$f, c (rf, c ) + $ab(r a b)$ac{rac)®bc(rbc) 



(50) 



However, we are interested in the density expansion of the free energy. The third cluster integral like all cluster 
integrals is a function of density n. We formulate now some general rules of picking out the leading terms of a density 
expansion of a certain cluster integral. A graph made of m points and I solid lines (i.e., it is a contribution to the 
mth cluster integral) is proportional to n m J dr\ . . . df m g l . Here the occurence of alltogether I Debye functions g 
with different arguments is scetched only schematically. Substituting Si — K?i one finds that the above graph gives a 
density contribution to the free energy density (with the additional factor <k, see Eq.(p5|)) of the order 



n ™ K 3 K -3m K l 



3/2 + 1/2 -m/2 



It should be mentioned that this simple estimate is valid only as far as the remaining dimensionless integral is 
convergent. Especially, it is not correct when bonds g p (r) of one argument r with p > 3 are present. These bonds 
are nonintegrable due to the singularity at short distances r. However, these short range divergencies are spurious. 
They are just an artefact of the expansion of Boltzmann factors in the Mayer-like graphs in powers of the screened or 
unscreened Coulomb potential. Since all Boltzmann factors are integrable at short distances it is always possible to 
eliminate these unphysical divergencies at any order by a suitable collection of the dangerous contributions EELEm , 
For example the g p (r) bonds (p > 3) may be resummed into the bond W — g 2 /2. At low densities, this bond can 
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be replaced by the bare truncated bond f(r) (see Eg. (p5|) ) . This is an integrable bond except for the logarithmic 
divergency at large distances. The latter divergency can be finally removed using the bond 



Ulab 

6 \ r 



-3x,r 



l)+/»6(r) , 



which leads to the appearance of the logarithmic term In k. This is just the content of the renormalization procedure 
employed in Eq.(J44|). It will be used in the remaining of this section. Nevertheless, the above estimation of the order 
of a cluster integral is correct as far as the leading contributions of the mth (to > 3) cluster integral are looked for. 
These contributions are obtained with bonds g and g 2 only for which the scaling analysis is applicable. One also proofs 
that the minimal number l m in{fn) of solid lines of a graph contributing to the mth cluster integral is 2to — [to/2] with 
[A] being the integer number of A. Therefore the leading contribution of the mth cluster integral is of the order of 



- ( :■! [ /;/ 



The leading contribution from the third cluster integral is the graph with 5 solid lines (stemming from the TL3 H3 H.4 
term, see Eq.(33)), 



S3 = ^3 
3 3! 




"o" X! c a c >> c c dr a df b df c g a b{r a b) -z 



a.b.c 



Jac(r a c)} 2 2 [9bc(rbc)} 2 



and gives a contribution to the order n 5 ' 2 



u Vn3 ^ , ,2 , 

»3 > J CaCbCclablaJ 



a,b,c 



(51) 



The numerical constant 63 (an integral representation was already given by Haga p6[) is calculated in Appendix 0. 
It is expressable through the Euler's dilogarithm Lia(x) = Y^=x x k /k 2 , \x\ < 1 and equals in accordance with 0] 



63 = 2tt^ 



12 



Li 2 - 



= 10.13477. 



(52) 



Following the general scheme the n 3 terms should be obtained from all graphs with 6 solid lines, i.e., from the diagrams 




and 




(53) 



representing the H4 H.4 TC4 and 7i 3 7i 4 7i 5 terms. 

However, the corresponding integrals are divergent. The divergencies occur at short distances and are connected 
with a singular accumulation of g 2 bonds or with the occurence of a bond g 3 . These short range divergencies 
are unphysical and can be removed by a renormalization procedure as discussed above. Regard therefore the two 
corrections to the leading contribution S3, 



n 3 ^ r 

SlSs = -ft - /, CaCbC c df a df b df c $ab(rab)$ac(r ac )®bc(rbc) ■ 
a,fe,c 



(54) 



and 



62S: 



2<->3 



n 

~2 



^c a c b c c dr a dr b dr c g ab (r ab ) 2 

„ t „ J 



^ac\ r ac) ' ydacvac) 



<f>bc(rbc) 



(55) 



and expand the corresponding contributions in terms of k. 

First, consider the term S1S3. As it is shown in Appendix H we get up to the order n 3 , 



S X S 3 = - Vn 3 Y^ c a c b c c l 2 ab l 2 ac l 2 c 



12 



• In kI„ 



C3c 



(56) 
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where the constant 



('3o 



71" 
12 



lim 167r 

(5^0 



l ' 5 dk 



arctan 



27r 4 hi(5 + C E 



2.89476 . 



has been introduced. Consider now the term 82S3. Up to the order n 3 we get (see Appendix ||) 



with 



5 2 S 3 = Vr? 2J CaCbCjabllj 
a,b,c 



2tt 2 



17 1 



— In nl ac + — I C E - yr + - In 3 ) In k1 



C3/3 



2tt 2 



C3/3 



121 
~L8~ 



C* s -l + ln3 
17 



C E -^ + ln3 
6 



In 3 [81n2-31n3]-81n 2 2 



2C% 



Notice that the leading term of 52S3 is of the order n In n. 



-C E _4Li 2 (4)-2Li 3 (-| 

3l„2, 



= 8.85348. 



(57) 



(58) 



(59) 



C. Fourth cluster integral 



The diagrammatic representation of the fourth cluster integral is found from the above diagrammatic rules and 
reads in accordance with 101 , 



54 = h 6 



12 



24 





+ 




+ 


1 


X 


1 


+ 12 


1 




r 


•c 












+ 12 




xy 


+ ... 


+ 3 






; 


'F == ^ 


1 




+ 6 


.( 


P^ 


+ ... 


■ 




c_ 




J 



(60) 



The dots beyond a diagram of a given type represent all diagrams constructed from the basic one by adding an 
arbitrary number of solid lines between points connected by 2 lines (for the case of outer lines) or by one line (for the 
case of inner lines). Introducing a generalized Mayer bond of second type, 



^ab(rab) = exp[- g a b(r a b)] - 1 



with the graphic representation, 



(61) 



one writes the fourth cluster integral in the following manner Pq , H ,B| , 



Si 




(62) 
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The leading contribution from the fourth cluster integral are the two graphs with 6 solid lines (stemming from the 
U-iU-iU-iU-i term), 



S 



4! 








(63) 



According to the rules formulated in the subsection pll F3| they are of the order of n 5 / 2 , and give [ p6|j25| 

4 

S\ = b 4 ^P e Yl c a c b c c c d e 3 a ele 3 c e 3 d , (64) 



with the numerical constant & 4 = b. 



4 — u 4,q T t>4,0i 



a,b,c,d 



— x l \ p — &12 I o — ^23 



h,a = 777 /dxidx 2 dx3 ( — 

16 J III y X\2 \ £23 



•i'3 



1 



-Xi „—X2 „ — X 3 „-X 12 f, — X 13 p — £23 



h.p = 777 /rfx 1( ix2rfx3 

24 7 xi £ 2 x 3 Xvx x i3 x 23 



The constant &4 >Q , is calculated in Appendix |Aj. It equals 



6 4q = 27T d 



7T 3 -. / 1\ 43 
T2 +Ll2 "3 + 3 ln 4 



8.05281 



(65) 
(66) 

(67) 



The constant b^p was evaluated numerically by a Monte Carlo calculation and equals, in closed correspondence with 



°4,/3 



1.785 ±0.001 . 



(68) 



Consider now all diagrams with 7 solid lines representing the H3 H.3 H4H4 and H3 H3 H3H5 terms and contributing 
to the order n 3 . 



SS 4 = I 12 




12 



O 



12 




12 




+ o(n 3 ) 



(69) 



Again the first diagram in the above equation is divergent and we may apply the renormalization procedure as 
described in the section III Bj. The details are given in Appendix O We obtain 



SSi = —f3 7 ^ c a c b c c c d < e a ( 

r, h ^ A V 



4 p 3 3 
b e c e d 



4tt 3 

-— In Kl ab - c 4 
y 



la&y 



d c 4e ) j 



a,b,c,d 

where the constants C4 = c 4a + C4/3 + C4 7 + C45 and C4 e has been introduced, 
"17 



Cic 



4tt 3 
~9 



6 



2CV 



16tt 2 '°° 



?r arctan — < 3 arctan — 

3 Jo (fc 2 + l) 2 2 \ 3 



In (k 2 + 9) 



c 4,/3 = o /rfx 1 dx 2 dx 3 



:Cl 



1 /" /e^ 1 

C4, 7 = 777 dxidx 2 dx 3 

16 J \ xi 



-X 2 p~X 3 p — Xi 2 p — Xi3 p -X 2 3 



X 2 X3 X\2 #13 X 2 3 



X3 



X\2 I X13 



c 4,s = 77 /dx!dx 2 (ix 3 

1 



X l / P~ X 3 



,—X\2 



2 _ 



^23 p— *13 



C4,. 



dxi(ix 2 (ix 3 



*i \ 13 / \ »12 / ^23 £13 

„— Xl \ 2 / e ~X 3 \ 2 e ~ x 12 g- a; 23 g-^2 



Xl 



X3 



X12 X23 X 2 



(70) 

(71) 
(72) 
(73) 
(74) 
(75) 
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The above multiple integrals for c^p, c^s and C4 i£ can be reduced to single or double integrals and are evaluated by 
numerical integration. The remaining integral is estimated by a Monte Carlo run. As a result we obtain: 

c 4 = 53.66 ±0.02, 
c 4 , e = 26.719... . (76) 

The fifth and sixth cluster integral contributing to the order n 3 are considered in the next section. 

D. Fifth and sixth cluster integral (leading terms) 

The diagrammatic representation of the fifth and sixth cluster integral is given by a huge number of basic diagrams. 
Thus the fifth cluster integral is represented in Rcf . || . According to our aim of presenting a virial expansion of the 
free energy it suffices to give the leading contributions to both cluster integrals. 

The leading contribution of the fifth cluster integral is given by the four graphs with 8 solid lines (stemming from 
the H4 {Ti-z) term), 



5 



h\ ls 




60 




-Q 



120 




(77) 



According to the rules formulated in the subsection IIIB they are of the order of n , and give 

,5 



S l 5 = c 5 -^P 8 Yl c a c b c c c d c f e a ele 3 c e 3 d e 3 f , 

a,b : c : d,f 



with the numerical constant C5 = c 5 a + C5 3 + c 5 7 + c 5 g. 

The four constants introduced here are determined by the dimensionless integrals. 



(78) 



1 / , 

C5 a = — /(ixidX2dx3(iX4 



x\ „—x 12 „-x 2 p — X 2 3 (,—X3 p—X3A p — Xi p—xa 



c 5,/3 



C5,7 



<ixi<ix2<ix3dx4 



X\ X\1 X 2 X 2 3 X3 X34 X4 Xil 
2 



-Xi 



J'l 



-X2 p — Xi p -X 13 p — X 3 4 p — X4,2 p~X23 



X 2 Xi X13 X34 X42 X 2 Z 



1 f ( e~ Xi 

— /(ixi(ix2dx3dx4 

16 J \ xi 



,-X 2 



-Xl3 p — X3 p — X3i / p — X24 



1 f fe- Xi 

C5,8 = ~. dx.1dx2dx.3dx4 

4J \ xi 



2 _ 
e 



X2 J X\Z X-i X 34 

X2 p-X\3 P ~X3 p—X34 



X2A 

2 



X2 X\S X 3 X 3 4 \ X 2 4 

A numerical evaluation (the "irreducible" integrals were estimated by Monte Carlo simulations) gives: 

c 5 = 168.2 ±0.6 . 
In the sixth cluster integral we can pick out the contributions stemming from the (TI3) term, 



(79) 
(80) 
(81) 
(82) 

(83) 



66 ~ 6! 



120 




180 




180 




10 




180 




(84) 



The five contributions are again of the order n 3 , and give 

,6 



4 = - c 6^/3 9 J2 c aCbC c c d c fCg e 3 a e 3 b e 3 c e 3 d e 3 f e 3 g 

a,b,c,d,f,g 



with the numerical constant cq = CQ. a 4- C6,/3 + c$ n + Cq^s + C6. e given by the five dimensionless integrals, 



(85) 
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1 f ( e~ Xl 

ce.a = 777 dxidx.2dx 3 dx4,dx 5 

48 J V xx 



1 f fe~ Xl 

cg,p = 777 dx 1 dx 2 dx 3 dx4 i dx 5 

16 J \ x 1 



£12 /p-£23 



-2=34 



-£45 



-x 5 



2^12 V ^23 / ^34 V 2: 45 / X 5 

£12 p~£23 / p-X34 \ 2 p-^45 p-£5 p-£52 



2 -„,,.. 

e ■"" e 



a?i2 2:23 \ ^34 / 2:45 x 5 x 52 



(86) 

(87) 



c 6-7 = o /rfxiiix2(ix3(ix4(ix5 
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-£12 p — £23 p-^34 p — £45 p-%5 n-XH3 p — £24 
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If e " L e ■"" e 
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C6,. 



-X\ p — £12 p — £23 p — £34 — £45 -15 — £51 — £24 „ — £3 



1 /" e ""■ e 

— /(ixi(ix2(ix3dx4(ix5 

4 7 *i *12 2^23 2:34 x 45 2:5 £51 x 24 2! 3 



(89) 
(90) 



The numeric evaluation results in: 



c 6 = 339.8 ±0.3. 



(91) 



According to the rules formulated in Sec. IIIB there are no other terms contributing to the order n . Thus in this 
section we have explicitely calculated all terms of the free energy of a classical plasma consisting of charges with the 
same sign and immersed in a neutralizing background up to the third order in density. The final result is given in 
Sec. fvl 



IV. FINAL RESULTS AND CONCLUSIONS 



The final expression for the interaction part of the Helmholtz free-energy density to the order n 3 is obtained 
from Eg. (35) by summing the corresponding contributions £2 (E q.jiq )), S3 (Eqs. (pi]) , ( 52 ) , (pq) , (p7[) , (pq) , (p9|) ) , S4 
(Eqs.@^),(f||,@,@), S 5 (Eqs.@,@), and S 6 (Eqs.Q,^) with the following result 



(3F C 
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-rf- = -777- _ n2 P 3 ^2 c a c b e 3 a el -lnnl ab + a 2 \ - n 2 K/3 4 y^ c a c b e a 
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a . b 
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b e c e d 



P / 4 c a c b c c e a e b e c 



a,b,c 



2ir 2 , 2 , 8tt 2 / 17 1 , . . , , 
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V/3 5 Y, c ac b el 



a . b 



— In Kl ab + c 2 



P / J c a c b c c e a e b e c 



1 k 2 (3 7 J2 c aCbC c c d \ e\444 



a.b.c.d 



4tt 3 



ItikL 



a,b,c 



(■i 



— \nnl ab + c 3 , a 



4444 C4e } 



C5 ~^ 8 5Z c a c b c c CdC}e A a e 3 b e 3 c e 3 d e 3 + c 6 —/3 9 V c a c b c c c d CfC g e 3 a ele 3 c e 3 d e 3 e 3 g + o (n 3 ) , (92) 

a,b,c,d,f a,b,c,d,f,g 



All constants in Eq. (|92|) are given in Sec. III. In the case of a one component plasma (OCP) consisting of one sort 
of ions with density n , charges Ze and moving in a neutralizing background all thermodynamical functions may be 
expressed in terms of the plasma parameter T = /3Z 2 e 2 (47rn/3) 1 ' 3 only. For the excess internal energy defined by 



u = (3 



d((3F c /N) 
dp 



we obtain from Eq.(|92|) 

u(r) = Po r 3/2 + Pl r 3 in r + p 2 r 3 + P3 r 9/2 in r + P4 r 9/2 + p 5 r 6 in 2 r + P6 r 6 in r + P7 r 6 , (93) 
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with the constants 



Po = 


~2 Vd ^ 




V2 = 


Q 


> 


VA = 


0.2350 




P6 = 


-2.0959 


) 



Pi 



>1 


= -r^ 




16 




81 


'b 


16 ' 


'7 


= 0.0676 



This expansion completes the result of Ref. pTj for the small-r expansion of the OCP excess internal energy. 

In conclusion, in this paper we have studied the low density expansion of the Helmholtz free energy of a classical 
system of pointlike ions embedded in a neutralizing background and interacting via Coulomb forces. Such a purely 
Coulomb description is appropriate for a large variety of physical systems. To obtain the virial expansion of a classical 
Coulomb system the Hubbard-Schofield transformation for the configurational integral in collective variables is used. 
The original Hamiltonian was maped onto an Ising-like Hamiltonian with harmonic and anharmonic contributions. 
The coefficients of the new effective Hamiltonian are given in terms of the structure factors of a reference system. In 
our case with pointlikc ions the corresponding reference system is an ideal gas system. 

However, it is possible to generalize the present approach and to consider mixtures of ions with internal structure, 
i.e., systems with additional short range interactions. In this case the corresponding reference system is a system 
with particles interacting solely via short range forces. By introducing short range repulsive interactions u ab (r) (for 
example hard core repulsions) it is also possible to describe systems of charges with different signs within the present 
approach. In this case the corresponding cluster expansion is built with the integrable bonds g a b(r), 9ab( r ) an d 
exp(-(lu ab (r)) exp(-g ab (r)) - 1 + g ab (r) - (l/2)g 2 ab (r) . 

Expanding the anharmonic contributions of the effective Hamiltonian to the configurational integral we have ob- 
tained the cluster expansion of a classical Coulomb gas (Eq.(p4[)). The exact density expansion of the free energy 
has been performed up to the order n 3 (Eq.(|92|)). The virial expansion of other thermodynamic functions can be 
obtained from (|92| ) using thermodynamic identities. Besides its own conceptual interest, the virial expansion might 
be useful in studying regimes with sufficiently low density and/or sufficiently high temperature. However, using virial 
expansions for practical calculations one has to take into account that a virial expansion is an asymptotic expansion 
of the free energy function. The convergence radius of this asymptotic expansion is unknown. Therefore it does not 
make sense to extent the region of the low density limit to higher densities by calculating higher order terms. We 
conclude therefore that the obtained virial expansion is applicable for practical calculations only in the low density 
region (or in the weak coupling regime). 

Finally, the present approach can be generalized to quantum plasmas. Thus one may employ Morita's method of 
effective potentials to map the quantum plasma system with bare Coulomb interaction to a classical system interacting 
via effective potentials. The latter split into a Coulomb part and a short range part. Further one expresses the 
Coulomb interactions through collective variables and maps the original quantum system via the Hubbard-Schofield 
transformation to a classical reference system with sole short range interactions. Alternatively one may also use the 
Feynman-Kac formalism. However, explicite calculations to the order n are much more complicated as in the classical 
case. 
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APPENDIX A: CALCULATION OF CONSTANT B 3 AND B 4 , a 
Consider the leading contribution to the third cluster integral in k space representation 
I 3 /T\ _ " 3 V(-1) 5 ^ _ , , 2 , 2 fdkidk 2 dk 3 4tt 



x^ c ChC i a 2 i 2 f dkldk: 

, J (2ir 
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3! // \A 2 2! 2! ^ — ~<- """" <">'-" J (2) 9 tf + K 
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Substituting Xj = fcj/« one expresses the leading term of the third cluster integral through the constant a\ with 

1 f dx\x\ 
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m) 
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Performing the integration over cc one obtains 



2tt 2 y 
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Now in Eq.(A2) the integration over x\ can be carried out and one obtains to] 



2n 2 



12 2 4 



1^(4) 



(A2) 



(A3) 



(A4) 



(A5) 



Using the properties of the Euler's dilogarithm one arrives at Eq.(|55_ 

Consider now the constant bj, a defined in Eq. (pq) . Using the Fourier representation of exp[— x]/x we get a similar 
representation as for 63 (Eq.(A2)), 



dx\x\ 



(A + iy 



I 2 (ti) 



with the same I{y) given by Eqs.(|A3|) and (A4). The remaining integration in Eq.(A6) leads to the result, 
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Further transformations lead to the more compact expression Eq.(|67|) 



(A6) 



(A7) 



APPENDIX B: CALCULATION OF 5iS 3 AND 82S3 

After integration over one volume integral and up to the order n 3 one represents Eq.(p3) in the following manner, 

Vn 3 



hS 3 = -g- Yl CaCbCc d ridr 2 - gl b {r 12 ) -gl c {ri) $ QC ( 



With the dimcnsionlcss variables Xi — kTj one has, 
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is representable as 
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J(x 2 ) 
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Within the order o{n 3 ) we can substitute the contributions with exp[— 3x 2 ]/3!2;2 . . . etc. by their small x 2 expansions 
to get a convergent expression, 
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We get with J(x 2 — > 0) = tt 3 /x 2 the following expression 
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■\2 



exp 



vacK 

X 2 



- 1 



x 2 2 x\ 



+ o(n 3 ) 



(B6) 



After integration and performing the limiting procedures we obtain the result Eq. ( p6| ) . 

Up to the order n 3 one represents Eq.(pq) after integration over one volume integral as follows, 



8 2 S 3 = ~ Vn 3 ^ c a CbC c dfidr 2 g a b(ri 2 ) - 5 2 c ( r i) 

a,b,c 



<£ac(r 2 ) - ^9 2 ac (r 2 ) 



+ o(n 3 ) 



(B7) 



Introducing dimcnsionless variables x% = nri and performing the integration over x\ one obtains, 

5 2 S 3 = -47T 2 Vn 3 y^c a c b cJ ab l 3 ac l 2 bc / dx 2 x 2 { e~ x * [Ei(-x 2 ) + ln3] - e* 2 Ei(-3x 2 )} 

/ K p-4x 2 



a,b,c 

— 3^2 



-1 e 

6 x 3 . 



24 x\ 



± 



+ o(n 3 ) 



(B8) 



The higher order terms containing exp[— 4x 2 ]/xf ... can be substituted by their small x 2 expansion to cut the 
singularity at small x 2 . One gets, 

6 2 S 3 = - 4ir 2 Vn 3 Y, CaCbcJabl 3 ac l 2 bc / dx 2 x 2 { e~ x * [Ei(-x 2 ) + In 3] - e* 2 Ei(-3x 2 )} — — -g- 



a,b,c 



6 x% 



2x 



LacK J- \vacfc) *• 



21 4! :4 



5! 



± ... (C £ -l + ln3x 2 ) + 0(7^ 



(B9) 



u 2 v. ^2 

Summing up the singular at x 2 = contributions a finite expression is obtained, 

5 2 S 3 = -^ 2 Vn 3 Y c aCbC c labllf bc / dafe&a { e~* 2 [Ei(-x 2 ) + In 3] - e* 2 Ei(-3x 2 )} 



-1 e 



-3a;2 



6 x 3 



^ x 2 



(C* £ -l + ln3x 2 



exp 



.1:2 



l>acK (,'ac K J ('ac K j 



•('2 



2x 2 



6xi 



+ o(n 3 ) . (BIO) 



The remaining convergent integral can be calculated by introducing an upper and lower bound and then by performing 
the limiting procedure to 00 and 0, respectively. After some lengthy calculations one arrives at Eq.(p8[). 



APPENDIX C: CALCULATION OF <5iS 4 



The first diagram in Eq.(p9|) represents a divergent contribution. Consider therefore the convergent correction to 
the leading contribution of the first diagram in Eq.(62) 
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5iS 4 = - < (i 




Within the order o (n 3 ) one obtains 

-.4 



<5i5 4 = — y^ c a c b c c c d df a df b df c dr d 2 

a,b,c,d 

9ac(r ac ) g bd {r bd ) -g 2 cd {r cd ) + o(n 



®ab(ra,b) ~ ^9l b (rab) 



or 



4^4 



— 3x 



<5i5 4 = - "T^x/? 7 X! c a c b c c c d e^ebe^c. d , ^ -~ 



5\x 5 



a,b,c,d 



± 



^lele d Jdx—J 4 {x) - 
. J J A (x^0) + o(n 3 ) 



>>ab K 



where 



is representable as 



Ji{xi) = / dx~2dx~3 — — 

J x 23 



— 22;23 g — X2 p — x 13 

x| 3 x 2 x 13 



J4(f) = T^y j d% 



(An) 2 Air k -z g 

,, 9 —. 9 — arctan - e 
(fc 2 + l) 2 k 2 



Taking into account that J4(0) = 87r 2 /3 we get the convergent expression 
5x54- 



5— 1 ) c a c b c c c d e a etele d \\m lim / dx — — J 4 (») - 

K z * — ' e^O R^oo J OX 

a.b.c.d 

1 + T ^ 2^~ + 6^3-} + »(»)■ 



8^ 2 



3(7a6«0 3 



exp 






After integration and performing the limiting procedures one gets the logarithmic contribution of Eq. 
constant C4 Q (Eq.(|7l|)). 



(CI) 



(C2) 



(C3) 



(C4) 



(C5) 



(C6) 
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